R Notebook

library(tidyverse)
library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

library(lubridate)

Attaching package: ‘lubridate’

The following object is masked from ‘package:cowplot’:

    stamp

The following objects are masked from ‘package:raster’:

    intersect, union

The following objects are masked from ‘package:dplyr’:

    intersect, setdiff, union

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

library(mgcv)

Loading required package: nlme

Attaching package: ‘nlme’

The following object is masked from ‘package:raster’:

    getData

The following object is masked from ‘package:dplyr’:

    collapse

This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

source("UVP_2017_library.R")

Loading required package: future

Attaching package: ‘future’

The following object is masked from ‘package:raster’:

    values

theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')

Particles Only

Read In Data

bes<- read_csv("dataOut/binned_EachSize.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.

bds <- read_csv("dataOut/binned_DepthSummary.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.

ues <- read_csv("dataOut/unbinned_EachSize.csv")

Parsed with column specification:
cols(
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = ""),
  depth = col_double(),
  psd_gam = col_double(),
  vol = col_double(),
  sizeclass = col_character(),
  lb = col_double(),
  ub = col_double(),
  binsize = col_double(),
  TotalParticles = col_double(),
  nparticles = col_double(),
  n_nparticles = col_double(),
  biovolume = col_double(),
  speed = col_double(),
  flux = col_double(),
  flux_fit = col_double(),
  GamPredictTP = col_double()
)

uds <- read_csv("dataOut/unbinned_DepthSummary.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.

PhoticBase <- 160
OMZBase <- 850

Summary Statistics

PlotNParticles <- uds %>% 
  ggplot(aes(x = tot_nparticles, y = depth, col = profile)) +
 facet_wrap(~project) +
 geom_point(alpha = 0.3, shape = 1) +
scale_y_reverse() + scale_x_log10()

PlotNParticles

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

#FSG4 <- gam(tot_nparticles~ s(depth, k = 3)  + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

summary(FSG1)

Family: gaussian 
Link function: identity 

Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, 
    bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.002      0.106   84.96   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value  
s(depth) 1.204  1.366 3.346  0.0874 .
s(Day)   1.512  1.761 3.001  0.0386 *
s(Hour)  1.579  2.000 2.713  0.0389 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.183   Deviance explained = 24.2%
GCV = 0.73885  Scale est. = 0.67365   n = 60

#summary(FSG2)
#summary(FSG3)
#summary(FSG4)

summary(FSG1)$r.sq - summary(FSG2)$r.sq

[1] 0.07590024

summary(FSG2)$r.sq - summary(FSG3)$r.sq

[1] 0.05267319

summary(FSG3)$r.sq

[1] 0.0542546

But there is between projects:

ProjGam <- gam(tot_nparticles~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))

summary(ProjGam)

Family: gaussian 
Link function: identity 

Formula:
tot_nparticles ~ s(depth, k = 3) + factor(project)

Parametric coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          8.9983     0.4445   20.24   <2e-16 ***
factor(project)P16  17.0001     1.3754   12.36   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
         edf Ref.df     F p-value  
s(depth)   1      1 3.552   0.064 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =    0.7   Deviance explained = 70.9%
GCV = 12.411  Scale est. = 11.855    n = 67

PlotPSDmany <- uds %>% 
  filter(project == "ETNP") %>%
  ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
 
  #geom_path(aes(x = psd_gam)) + 
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
  geom_point(alpha = .6, size = 2, stroke = 1) +
  scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
  scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  labs(x = "Depth (m)", y = "Particle Size Distribution Slope")

PlotParticlesmany <- uds %>% 
  filter(project == "ETNP") %>%
  ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
 
  #geom_path(aes(x = psd_gam)) + 
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
  geom_point(alpha = .6, size = 2, stroke = 1) +
  scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
  scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  scale_x_log10() + theme(legend.position = "none") +
  labs(x = "Depth (m)", y = "Particles / L")

PlotFluxmany <- uds %>% 
  filter(project == "ETNP") %>%
  ggplot(aes(x = tot_flux_fit, y = depth, shape = factor(day(time)), fill = hour(time))) +
 
  #geom_path(aes(x = psd_gam)) + 
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
  geom_point(alpha = .6, size = 2, stroke = 1) +
  scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
  scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  scale_x_log10() + theme(legend.position = "none")



plot_grid(
  PlotParticlesmany,
  PlotPSDmany,
  rel_widths = c(2, 3)
  )

ggsave("figures/ParticlesPSDMany.png")

Saving 12 x 7.41 in image

ggsave("figures/ParticlesPSDMany.svg")

Saving 12 x 7.41 in image

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG4 <- gam(psd~ s(depth, k = 3)  + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

summary(FSG1)

Family: gaussian 
Link function: identity 

Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.96083    0.01988  -199.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F p-value    
s(depth) 1.713  1.917 56.197 5.5e-15 ***
s(Day)   1.000  1.000  1.078  0.3036    
s(Hour)  1.604  2.000  3.049  0.0285 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.635   Deviance explained = 66.1%
GCV = 0.026006  Scale est. = 0.023702  n = 60

#summary(FSG2)
#summary(FSG3)
summary(FSG4)

Family: gaussian 
Link function: identity 

Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.96083    0.01995  -198.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.715  1.919 55.141 7.66e-15 ***
s(Hour)  1.565  2.000  2.497   0.0488 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.632   Deviance explained = 65.2%
GCV = 0.02572  Scale est. = 0.023885  n = 60

summary(FSG1)$r.sq - summary(FSG2)$r.sq

[1] 0.03804023

summary(FSG2)$r.sq - summary(FSG3)$r.sq

[1] -0.004905911

summary(FSG3)$r.sq

[1] 0.6015425

Not a significant difference in PSD with respect to time.

But there is between projects:

ProjGam <- gam(psd~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))

summary(ProjGam)

Family: gaussian 
Link function: identity 

Formula:
psd ~ s(depth, k = 3) + factor(project)

Parametric coefficients:
                   Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        -3.96167    0.02491 -159.029   <2e-16 ***
factor(project)P16 -0.19977    0.07708   -2.592   0.0118 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F  p-value    
s(depth) 1.224  1.398 34.54 2.26e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.437   Deviance explained = 45.6%
GCV = 0.039117  Scale est. = 0.037234  n = 67

I wonder if I can show that the profiles aren’t statistically significanlty different. Or that they are for that matter… I think in that case, I run a gam with and without a parameter for profile… And then quantify the effect size of that parameter

Or follow this Gavin Simpson Post https://fromthebottomoftheheap.net/2017/10/10/difference-splines-i/

or anova.gam {mgcv}

Calculate gams for each profile, and then run anova.gam to see if they are different…

PlotNParticlesEP <- uds %>% 
  filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
 geom_point(alpha = 0.7, size = 2, stroke = 1) +
  #geom_path(aes(x = tot_nparticles)) +
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(2500, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
  labs(x = "Particles/L", y = "Depth (m)") +
  theme(legend.position = "none") +
  scale_shape_manual(values = c(1:5)) 

PlotNParticlesEP

I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)

PlotPSDEP <- uds %>% 
  filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
 geom_point(alpha = 0.7, size = 2, stroke = 1) +
  geom_path(aes(x = psd_gam)) +
  geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(2500, 0)) + scale_color_manual(values = c("gray20", "brown"))  +
  scale_shape_manual(values = c(1:5)) + labs(y = "", x = "Particle Size Distribution Slope")

PlotPSDEP

I may just cow these togther.

plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))

Removed 611 rows containing missing values (geom_point).Removed 611 rows containing missing values (geom_point).Removed 611 row(s) containing missing values (geom_path).

ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")

Saving 10 x 4 in image

ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")

Saving 10 x 4 in image

mainParticleComponents <- bds %>%
  filter(profile %in% c("stn_043", "p16n_100")) %>%
  select(project, profile, depth,
         tot_nparticles, small_nparticles, big_nparticles,
         tot_psd = psd, small_psd, big_psd,
         tot_flux_fit, small_flux_fit, big_flux_fit) %>%
  pivot_longer(cols = -c("project", "profile", "depth")) %>%
  separate(name, c("size", "meas")) %>%
  mutate(meas = recode(meas, nparticles = "particles/L")) %>%
  mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘select’ for signature ‘"spec_tbl_df"’

Flux small and flux tot track so closely because ag > psd. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top.

eg_dataline <- bds %>% 
  filter(profile == "stn_043", depth == 162.5)
eg_slope =  eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)

eg_datablock <- bes %>%
  filter(profile == "stn_043", depth == 162.5)


eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)

eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)


EgNNP <- eg_datablock %>%
  ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() + 
  geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)", x = "Size (mm)")

EgNP <- eg_datablock %>%
  ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() + 
  geom_path(data = eg_df) + labs(y = "Normalized Particles" , x = "Size (mm)")

EgTP <- eg_datablock %>%
  ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() + 
  geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)", x = "Size (mm)")

plot_grid(EgNNP, EgTP, labels = c("A", "B"))

Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis

ggsave("figures/ExamplePSD163m.png")

Saving 7.29 x 4.5 in image

ggsave("figures/ExamplePSD163m.svg")

Saving 7.29 x 4.5 in image

Smooth flux and um disaggregation.

bds %>% 
  ggplot(aes(y = depth, x = Flux_Smooth, col = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10()

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

FSG1 <- gam(Flux_Smooth~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG2 <- gam(Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

FSG3 <- gam(Flux_Smooth ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

summary(FSG1)

Family: gaussian 
Link function: identity 

Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, 
    bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  34.7319     0.9174   37.86   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.878  1.985 12.389 8.56e-05 ***
s(Day)   1.881  1.985  3.235   0.0416 *  
s(Hour)  1.454  2.000  2.926   0.0233 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.352   Deviance explained =   41%
GCV =  56.33  Scale est. = 50.497    n = 60

summary(FSG2)

Family: gaussian 
Link function: identity 

Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  34.7319     0.9663   35.94   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.867  1.982 11.055 0.000216 ***
s(Day)   1.814  1.965  2.125 0.128023    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.282   Deviance explained = 32.6%
GCV =  60.77  Scale est. = 56.029    n = 60

summary(FSG3)

Family: gaussian 
Link function: identity 

Formula:
Flux_Smooth ~ s(depth, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  34.7319     0.9944   34.93   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F  p-value    
s(depth) 1.857  1.979 10.73 0.000277 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.239   Deviance explained = 26.3%
GCV = 62.291  Scale est. = 59.325    n = 60

summary(FSG1)$r.sq - summary(FSG2)$r.sq

[1] 0.0709392

summary(FSG2)$r.sq - summary(FSG3)$r.sq

[1] 0.04226961

summary(FSG3)$r.sq

[1] 0.2392424

bds %>% filter(project == "ETNP") %>% select(profile, depth, Flux_Smooth) %>% pivot_wider(names_from = profile, values_from = Flux_Smooth)

Something is off. All of the flux profiles are identical. Skip this

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
plt1 <- bds %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = DFP, col = factor(time), shape = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + xlim(c(0.5, 1.5))+ geom_vline(xintercept = 1) +
   scale_color_manual(values = c(rep("black", 5), rep("blue", 5))) + scale_shape_manual(values = rep(1:5, 2))

plotly::ggplotly(plt1)

What the heck is going on with DFP here. Why is it usually > 1 shouldn’t it be less than 1 when flux is decreasing? This very deep increasing flux seems improbable to me. Lets check the smooths. Or only go to 1000m.

, legend.background = element_blank(), legend.box.background = element_rect()

scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time)))  + geom_point(size = 3, stroke = 1)+
  geom_path() +
  scale_y_reverse(limits = c(1000, 0))+
  scale_x_log10(label = scientific_10) +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + 
  scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
  geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE) +
  theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm"))



pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlxLegend <- get_legend(pltFlx)

Removed 14 rows containing missing values (geom_point).Removed 14 row(s) containing missing values (geom_path).

pltFlx

#plotly::ggplotly(plt1)

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 3, stroke = 1)+
  geom_path() +
  scale_y_reverse()+
  #scale_x_log10() +
  scale_x_log10(breaks = c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20)), limits = c(20, 180)) +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) + 
  scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux", y = "Depth") + theme(legend.position = "none")

pltFlxZoom

#plotly::ggplotly(plt1)

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time)))  + geom_point(size = 3, stroke = 1)+
  geom_path() +
  scale_y_reverse(limits = c(1000, 0))+
  scale_x_continuous(limits = c(-2.1, .6), breaks = seq(from = -2, to = .75, by = 0.5)) +
  #scale_x_log10() +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + 
  scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  geom_vline(xintercept = 0) +
  labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")
  #labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")

pltDelta3

#plotly::ggplotly(plt1pos)

# #plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
# 
# pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
# 
# pgTop <- plot_grid(pltFlxNoLegend, pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
# pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""), label_size = 14)
# pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
# 
# pgBoth
# 
# ggsave("figures/FluxDeepDive.png")
# ggsave("figures/FluxDeepDive.svg")

Within panel drawing

pgTop <- ggdraw(pltFlxNoLegend 
       ) +
  draw_plot(pltFlxZoom, .4, .25, .55, .60) +
  draw_plot_label(
    c("","B"),
    c(.05, 0.55),
    c(1, 0.85),
    size = 16
  )

Removed 14 rows containing missing values (geom_point).Removed 14 row(s) containing missing values (geom_path).

pgTop

pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c(“C”, ""), label_size = 14)

I don’t know whats going on below here

pgBottom <- pltDelta3  + geom_rect(aes(xmin = -2, xmax = -1.15, ymin = 170, ymax = 1000), colour = "gray50", fill = NA, inherit.aes = FALSE) + draw_plot(pltFlxLegend , -1.9, -575, .7)
pgBoth <- plot_grid(pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
                    pgBottom + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
                    ncol = 1, rel_heights = c(4, 4), labels = c("A", "C"), label_size = 16)

Removed 33 rows containing missing values (geom_point).Removed 33 row(s) containing missing values (geom_path).

pgBoth

ggsave("figures/FluxDeepDive.png")

Saving 5 x 9 in image

ggsave("figures/FluxDeepDive.svg")

Saving 5 x 9 in image

# #plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
# 
# pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
# 
# pgTop <- plot_grid(pltFlxNoLegend + ylim(c(1000, 0)), pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
# pgBottom <- plot_grid(pltDelta3 + ylim(c(1000, 0)), pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""))
# pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
# 
# pgBoth
# 
# #ggsave("figures/FluxShallowDive.png")
# #ggsave("figures/FluxShallowDive.svg")

Test for day to day and hourly variability

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

DFG1 <- gam(DF/DZ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

DFG2 <- gam(DF/DZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

DFG3 <- gam(DF/DZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

summary(DFG1)

Family: gaussian 
Link function: identity 

Formula:
DF/DZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.1001     0.0173  -5.785  3.8e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.907  1.991 22.926 8.57e-08 ***
s(Day)   1.908  1.991  6.576  0.00212 ** 
s(Hour)  1.329  2.000  2.351  0.03593 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.502   Deviance explained = 54.5%
GCV = 0.020011  Scale est. = 0.017961  n = 60

summary(DFG2)

Family: gaussian 
Link function: identity 

Formula:
DF/DZ ~ s(depth, k = 3) + s(Day, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.10008    0.01804  -5.547  8.5e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.897  1.989 21.470 2.11e-07 ***
s(Day)   1.903  1.991  4.847  0.00988 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.458   Deviance explained = 49.3%
GCV = 0.021234  Scale est. = 0.019535  n = 60

summary(DFG3)

Family: gaussian 
Link function: identity 

Formula:
DF/DZ ~ s(depth, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.10008    0.01944  -5.149 3.37e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F  p-value    
s(depth) 1.878  1.985 19.24 9.12e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.371   Deviance explained = 39.1%
GCV = 0.023812  Scale est. = 0.02267   n = 60

summary(DFG1)$r.sq - summary(DFG2)$r.sq

[1] 0.04364864

summary(DFG2)$r.sq - summary(DFG3)$r.sq

[1] 0.08693722

summary(DFG3)$r.sq

[1] 0.3713016

png(filename = “./figures/CombinedP2Info.png”, width = 10, height = 8, units = “in”, res = 200) StationInfoPlot() dev.off()

#plot.new()
FluxGamPlot <- function(){
  par(mfrow = c(2,2))
  plot(DFG1)
  mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(1,1))
  mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(1,2))
  mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
}

FluxGamPlot()

png(filename = "./figures/FluxGamPlot.png", width = 10, height = 8, units = "in", res = 200)

FluxGamPlot()
dev.off()

png 
  2 

#Osps

(u mol C / m^3 / day)

bds %>% filter(project == "ETNP") %>%
  ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
  scale_x_continuous(limits = c(-1, 1)) +
  geom_vline(xintercept = 0) +   scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = bquote("Observed - Modeled Small Particle Flux"~(μmol/m^3/day))) +
  scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue")

#ggsave("..figures/FluxSizeShift.svg"

 ggsave("figures/FluxSizeShift.png")

Saving 6 x 4 in image

 ggsave("figures/FluxSizeShift.svg")

Saving 6 x 4 in image

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))

summary(OZG1)

Family: gaussian 
Link function: identity 

Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.024817   0.002255   11.01 1.77e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
            edf Ref.df     F  p-value    
s(depth) 1.8900  1.988 8.511 0.000405 ***
s(Day)   1.9131  1.992 6.058 0.003482 ** 
s(Hour)  0.5495  2.000 0.407 0.235090    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.341   Deviance explained =   39%
GCV = 0.00033503  Scale est. = 0.00030514  n = 60

summary(OZG2)

Family: gaussian 
Link function: identity 

Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.024817   0.002272   10.92 2.06e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F  p-value    
s(depth) 1.889  1.988 8.417 0.000432 ***
s(Day)   1.910  1.992 5.671 0.004881 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.331   Deviance explained = 37.4%
GCV = 0.00033672  Scale est. = 0.00030979  n = 60

summary(OZG3)

Family: gaussian 
Link function: identity 

Formula:
ospsDZ ~ s(depth, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.024817   0.002479   10.01 3.55e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F  p-value    
s(depth) 1.872  1.984 7.338 0.000983 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.204   Deviance explained = 22.9%
GCV = 0.00038738  Scale est. = 0.00036884  n = 60

summary(OZG1)$r.sq - summary(OZG2)$r.sq

[1] 0.01003818

summary(OZG2)$r.sq - summary(OZG3)$r.sq

[1] 0.1274581

summary(OZG3)$r.sq

[1] 0.2038316

OSMSGamPlot <- function(){
  par(mfrow = c(1,2))
  plot(OZG2)
  mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(1,1))
  mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
}

OSMSGamPlot()

png(filename = "./figures/OSMSGamPlot.png", width = 10, height = 6, units = "in", res = 200)

OSMSGamPlot()

dev.off()

png 
  2 

bds %>% filter(project == "P16") %>%
  ggplot(aes(y = depth, x = ospsDZ)) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(500, 0)) + geom_vline(xintercept = 0)

Error in eval(lhs, parent, parent) : object 'bds' not found

Trap data

Example particle size distribution

horizontalGamPlot <- dataGamHorizontal %>% ggplot(aes(x = resp_fit, y = depth, col = log(lb), group = lb)) + scale_y_reverse() + geom_point() + scale_x_log10(limits = c(10^-8, NA)) + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + geom_errorbar(aes(xmin = resp_lower, xmax = resp_upper), width = 10, alpha = 0.5)+ theme_bw()

Weber Bianchi Figs

Average of everything

Just 043

Fig 5 WB

bds %>% filter(profile == “stn_043”, depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = “darkgreen”) + geom_hline(yintercept = 950, color = “darkblue”) + geom_ribbon(alpha = 0.2) + labs(x = “PSD slope”)

All of them

043 only

Small particles biomass

For tom and danielle

P16 Flux

EK60